home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Languguage OS 2
/
Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO
/
gnu
/
libg_261.zip
/
libg_261
/
libg++
/
etc
/
fib
/
fib.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1994-05-13
|
14KB
|
616 lines
#include <stream.h>
#include <Integer.h>
#include <builtin.h>
// Fun and games with the Fibonacci function.
// Also demonstrates use of the `named return value' extension.
// fib1 and fib2 from Doug Lea, others from Doug Schmidt.
/**********************************************************************/
// Standard iterative version:
// Schmidt uses convention that fib(0) = 1, not fib(0) = 0;
// so I just increment n here.
Integer
fib1(int n)
{
++n;
if (n <= 0)
return 0;
else
{
Integer f = 1;
Integer prev = 0;
while (n > 1)
{
Integer tmp = f;
f += prev;
prev = tmp;
--n;
}
return f;
}
}
/**********************************************************************/
// n
// via transformed matrix multiplication of ( 0 1 )
// ( 1 1 )
// simplified knowing that the matrix is always of
// the form (a b)
// (b c)
// (where (a, b) and (b, c) are conseq pairs of conseq fib's;
// further simplified by realizing that c = a+b, so c is only
// calculated implicitly.
// Given all this, do the power via the standard russian
// peasant divide-and-conquer algorithm, with
// the current multiplier held in (p q)
// (q -)
//
// This example also shows that using procedure calls instead of operators
// can be faster for Integers
// operator version
Integer
fib2(int n)
{
++n; // 1-based for compatability; see above
Integer a = 0;
if (n <= 0)
return a;
else
{
Integer b = 1;
Integer p = 0;
Integer q = 1;
for(;;)
{
if (n == 1)
return a * p + b * q;
else if (odd (n))
{
Integer aq = a * q;
Integer bq = b * q;
a = a * p + bq;
b = b * p + aq + bq;
}
Integer qq = q * q;
q = (q * p) * 2L + qq;
p = (p * p) + qq;
n >>= 1;
}
}
}
// with all expressions named, including return value
Integer
fib2a(int n)
{
Integer a = 0;
++n;
if (n <= 0)
return a;
else
{
Integer b = 1;
Integer p = 0;
Integer q = 1;
for(;;)
{
if (n == 1)
{
Integer bq = b * q;
a *= p ;
a += bq;
return a;
}
else if (odd (n))
{
Integer aq = a * q;
Integer bq = b * q;
a *= p;
a += bq;
b *= p;
b += bq;
b += aq;
}
Integer qq = q * q;
Integer pp = p * p;
Integer pq = p * q;
q = pq;
q += pq;
q += qq;
p = pp;
p += qq;
n >>= 1;
}
}
}
// procedure call version
Integer
fib2b(int n)
{
Integer a = 0;
++n;
if (n <= 0)
return a;
else
{
Integer b = 1;
Integer p = 0;
Integer q = 1;
for(;;)
{
if (n == 1)
{
Integer bq; mul(b, q, bq);
mul(a, p, a);
add(a, bq, a);
return a;
}
else if (odd (n))
{
Integer aq; mul(a, q, aq);
Integer bq; mul(b, q, bq);
mul(a, p, a);
mul(b, p, b);
add(a, bq, a);
add(b, bq, b);
add(b, aq, b);
}
Integer qq; mul(q, q, qq);
Integer pp; mul(p, p, pp);
Integer pq; mul(p, q, pq);
add(pq, pq, q);
add(q, qq, q);
add(pp, qq, p);
n >>= 1;
}
}
}
// procedure call version, reusing variables
Integer
fib2c(int n)
{
Integer a = 0;
++n;
if (n <= 0)
return a;
else
{
Integer b = 1;
Integer p = 0;
Integer q = 1;
Integer bq, aq, qq, pp, pq;
for(;;)
{
if (n == 1)
{
mul(b, q, bq);
mul(a, p, a);
add(a, bq, a);
return a;
}
else if (odd (n))
{
mul(a, q, aq);
mul(b, q, bq);
mul(a, p, a);
mul(b, p, b);
add(a, bq, a);
add(b, bq, b);
add(b, aq, b);
}
mul(q, q, qq);
mul(p, p, pp);
mul(p, q, pq);
add(pq, pq, q);
add(q, qq, q);
add(pp, qq, p);
n >>= 1;
}
}
}
/**********************************************************************/
// Ullman memoizers:
class Ullman_Array
{
// Ullman's array implementation allows Initialization, Store, and Fetch
// in O(1) time. Although it takes O(n) space the time to initialize enables
// us to compute a Fibonacci number in O(log n) time!
// The basic concept of Ullman's array implementation is the use of a
// Hand_Shake_Array and an Index_Array. An Index location in the array is
// only considered initialized if the contents in the Hand_Shake_Array and
// Index_Array point to each other, i.e. they ``shake hands!''
private:
int max_array_range;
int max_set_size;
int *index_array;
int *hand_shake_array;
Integer *storage_array;
int current_max;
public:
Integer uninitialized;
Ullman_Array (int array_range, int set_size);
void store (int index, const Integer &item);
Integer fetch (int index);
};
inline
Ullman_Array::Ullman_Array (int array_range, int set_size)
{
max_array_range = array_range;
max_set_size = set_size;
index_array = new int[max_array_range + 1];
hand_shake_array = new int[max_set_size];
storage_array = new Integer[max_set_size];
current_max = -1;
uninitialized = 0; // No fibonacci number has value of 0.
}
// Store Item at the proper Index of the Ullman array.
// The Hand_Shake_Array determines whether the Index has already
// been stored into. If it has NOT, then a new location for it
// is set up in the Storage_Array and the Hand_Shake_Array and Index_Array
// are set to point at each other.
inline void
Ullman_Array::store (int index, const Integer &item)
{
int hand_shake_index = index_array[index];
if (hand_shake_index > current_max || hand_shake_index < 0
|| index != hand_shake_array[hand_shake_index])
{
hand_shake_index = ++current_max;
hand_shake_array[hand_shake_index] = index;
index_array[index] = hand_shake_index;
}
storage_array[hand_shake_index] = item;
}
// Returns uninitialized if not initialized, else returns Item at Index.
inline Integer
Ullman_Array::fetch(int index)
{
int hand_shake_index = index_array[index];
if (hand_shake_index > current_max || hand_shake_index < 0
|| index != hand_shake_array[hand_shake_index])
return uninitialized;
else
return storage_array[hand_shake_index];
}
/**********************************************************************/
class Memo_Fib : public Ullman_Array
{
public:
Memo_Fib (int fib_num);
Integer fib3 (int n);
Integer fib3a (int n);
};
// The total number of Items computed by the Fib routine is bounded by
// 4 * ceiling of log base 2 of Fib_Num. Of course, the basis of the
// recurrence is Fib(0) == 1 and Fib(1) == 1!
Memo_Fib::Memo_Fib (int fib_num) : Ullman_Array(fib_num, (4 * lg (fib_num)))
{
store (0, 1);
store (1, 1);
}
// Uses the memoization technique to reduce the time complexity to calculate
// the nth Fibonacci number in O(log n) time. If the value of ``n'' is
// already in the table we return it in O(1) time. Otherwise, we use the
// super-nifty divide-and-conquer algorithm to break the problem up into
// 4 pieces of roughly size n/2 and solve them recursively. Although this
// looks like an O(n^2) recurrence the memoization reduces the number of
// recursive calls to O(log n)!
Integer
Memo_Fib::fib3 (int n)
{
Integer fib = fetch(n);
if (fib == uninitialized)
{
int m = n >> 1;
fib = fib3 (m) * fib3 (n - m) + fib3 (m - 1) * fib3 (n - m - 1);
store (n, fib);
}
return fib;
}
// The same, with procedure calls instead of operators
Integer
Memo_Fib::fib3a (int n)
{
Integer fib = fetch(n);
if (fib == uninitialized)
{
int m = n >> 1;
Integer tmp;
mul(fib3a(m), fib3a(n - m), fib);
mul(fib3a(m - 1), fib3a(n - m - 1), tmp);
add(fib, tmp, fib);
store (n, fib);
}
return fib;
}
/**********************************************************************/
// Here's a linear-time dynamic programming solution to the same problem.
// It makes use of G++/GCC dynamic arrays to build a table of ``n'' partial
// solutions. Naturally, there is an O(1) space solution, but this O(n)
// approach is somewhat more intuitive and follows the ``original'' recurrence
// more closely.
Integer
fib4 (int n)
{
#if defined (__GNUC__) && !defined (__STRICT_ANSI__)
Integer table[n + 1];
#else
Integer *table = new Integer[n + 1];
#endif
table[0] = 1;
table[1] = 1;
for (int i = 2; i <= n; i++)
table[i] = table[i - 1] + table[i - 2];
#ifdef __GNUC__
return table[n];
#else
Integer tmp = table[n];
delete [] table;
return tmp;
#endif
}
/**********************************************************************/
// The extended integers provide numbers of the form:
// (base+sqrt(5)*Rad_5)/2^Pow_2
// These are used to find the solution to the closed form of fib(n).
struct Extended_Int
{
Integer base;
Integer rad_5;
int pow_2;
Extended_Int (const Integer ¶m1, const Integer ¶m2, int param3);
Extended_Int (const Extended_Int& param);
Extended_Int (void) {}
friend inline Extended_Int operator- (const Extended_Int ¶m1,
const Extended_Int ¶m2);
friend Extended_Int operator* (const Extended_Int ¶m1,
const Extended_Int ¶m2);
friend inline Extended_Int sqrt (const Extended_Int ¶m1);
friend Extended_Int pow (const Extended_Int ¶m1, int num);
};
inline
Extended_Int::Extended_Int (const Integer ¶m1, const Integer ¶m2,
int param3)
{
base = param1;
rad_5 = param2;
pow_2 = param3;
}
inline
Extended_Int::Extended_Int (const Extended_Int ¶m)
{
base = param.base;
rad_5 = param.rad_5;
pow_2 = param.pow_2;
}
inline Extended_Int
operator- (const Extended_Int ¶m1, const Extended_Int ¶m2)
{
Extended_Int temp;
temp.base = param1.base - param2.base;
temp.rad_5 = param1.rad_5 - param2.rad_5;
temp.pow_2 = param1.pow_2;
return temp;
}
Extended_Int
operator* (const Extended_Int ¶m1, const Extended_Int ¶m2)
{
Extended_Int temp;
temp.base = param1.base * param2.base + 5L * param1.rad_5 * param2.rad_5;
temp.rad_5 = param1.base * param2.rad_5 + param1.rad_5 * param2.base;
temp.pow_2 = param1.pow_2 + param2.pow_2;
while (temp.pow_2 > 0 && !(odd (temp.base) || odd (temp.rad_5)))
{
temp.base >>= 1;
temp.rad_5 >>= 1;
temp.pow_2--;
}
return temp;
}
inline Extended_Int
sqrt (const Extended_Int ¶m1)
{
return param1 * param1;
}
Extended_Int
pow (const Extended_Int ¶m1, int num)
{
if (num > 1)
return odd (num)
? param1 * sqrt (pow (param1, num >> 1)) : sqrt (pow (param1, num >> 1));
else
return param1;
}
/**********************************************************************/
// Calculates fib (n) by solving the closed form of the recurrence.
class Closed_Form : private Extended_Int
{
private:
Extended_Int cons1;
Extended_Int cons2;
public:
Closed_Form (void): cons1 (1, 1, 1), cons2 (1, -1, 1) {}
Integer fib5 (int n);
};
Integer
Closed_Form::fib5 (int num)
{
Extended_Int temp = pow (cons1, num + 1) - pow (cons2, num + 1);
while (temp.pow_2 > 0 && !(odd (temp.base) || odd (temp.rad_5)))
{
temp.base >>= 1;
temp.rad_5 >>= 1;
temp.pow_2--;
}
return temp.rad_5;
}
/**********************************************************************/
static const int DEFAULT_SIZE = 10000;
void dofib1(int fib_num)
{
start_timer ();
Integer result = fib1 (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib1 = " << result << ". Time = " << time << "\n";
}
void dofib2(int fib_num)
{
start_timer ();
Integer result = fib2 (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib2 = " << result << ". Time = " << time << "\n";
}
void dofib2a(int fib_num)
{
start_timer ();
Integer result = fib2a(fib_num);
double time = return_elapsed_time(0.0);
cout << "fib2a = " << result << ". Time = " << time << "\n";
}
void dofib2b(int fib_num)
{
start_timer ();
Integer result = fib2b(fib_num);
double time = return_elapsed_time(0.0);
cout << "fib2b = " << result << ". Time = " << time << "\n";
}
void dofib2c(int fib_num)
{
start_timer ();
Integer result = fib2c(fib_num);
double time = return_elapsed_time(0.0);
cout << "fib2c = " << result << ". Time = " << time << "\n";
}
void dofib3(int fib_num)
{
start_timer ();
Memo_Fib Memo_Test (fib_num);
Integer result = Memo_Test.fib3 (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib3 = " << result << ". Time = " << time << "\n";
}
void dofib3a(int fib_num)
{
start_timer ();
Memo_Fib Memo_Test (fib_num);
Integer result = Memo_Test.fib3a (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib3a = " << result << ". Time = " << time << "\n";
}
void dofib4(int fib_num)
{
start_timer ();
Integer result = fib4 (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib4 = " << result << ". Time = " << time << "\n";
}
void dofib5(int fib_num)
{
Closed_Form Rec_Test;
Integer result = Rec_Test.fib5 (fib_num);
double time = return_elapsed_time(0.0);
cout << "fib5 = " << result << ". Time = " << time << "\n";
}
int
main (int argc, char *argv[])
{
int fib_num = (argc > 1) ? atoi (argv[1]) : DEFAULT_SIZE;
cout << "Results for fib(" << fib_num << "):\n\n";
dofib1(fib_num);
dofib2(fib_num);
dofib2a(fib_num);
dofib2b(fib_num);
dofib2c(fib_num);
dofib3(fib_num);
dofib3a(fib_num);
// dofib4(fib_num); // This uses too much mem for large n!
dofib5(fib_num);
return 0;
}